{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: Enter the code block and hit Ctrl+Enter to execute it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p: 3242.546045484618 Pa\n",
      "B_2: -2.238945068494697e-05 m^3/mol^1\n",
      "B_3: 1.0371257814774649e-09 m^6/mol^2\n",
      "B_4: 3.340005219645591e-14 m^9/mol^3\n",
      "betaV: 10.809571850439895 Pa/K\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import multicomplex as mcx\n",
    "from scipy.special import factorial\n",
    "\n",
    "R = 8.314462618 # J/mol/K\n",
    "\n",
    "# Values for argon\n",
    "Tc = 150.687 # K\n",
    "pc = 4863000.0 # Pa\n",
    "\n",
    "class vdWEOS:\n",
    "    a = (27/64)*(R*Tc)**2/pc\n",
    "    b = (1/8)*(R*Tc)/pc\n",
    "\n",
    "    def alphar(self, T, rho): \n",
    "        return -np.log(1.0-self.b*rho)-self.a*rho/(R*T)\n",
    "\n",
    "# Instance of the class\n",
    "vdW = vdWEOS()\n",
    "\n",
    "# This state point\n",
    "T = 300 # K\n",
    "rho = 1.3 # mol/m^3\n",
    "\n",
    "# Pressure\n",
    "dalphardrho = mcx.diff_mcx1(lambda r: vdW.alphar(T, r), rho, 1, True)[1]\n",
    "p = rho*R*T*(1+rho*dalphardrho)\n",
    "print(f'p: {p} Pa')\n",
    "\n",
    "# The first few virial coefficients\n",
    "# ...\n",
    "# Derivatives of alphar at zero density up to fourth order\n",
    "ders = mcx.diff_mcx1(lambda r: vdW.alphar(T, r), 0.0, 4, True)\n",
    "for n in range(2,5):\n",
    "    B_n = ders[n-1]/factorial(n-2)\n",
    "    units = f'm^{3*(n-1)}/mol^{(n-1)}'\n",
    "    print(f'B_{n}: {B_n} {units}')\n",
    "\n",
    "# Isochoric tension\n",
    "# ...\n",
    "# First derivative of alphar w.r.t. T at constant density\n",
    "dalphardT = mcx.diff_mcx1(lambda T_: vdW.alphar(T_, rho), T, 1, True)[1]\n",
    "# Cross (density, temperature) derivative of alphar\n",
    "d2alphardTdrho = mcx.diff_mcxN(lambda x: vdW.alphar(x[0], x[1]), [T, rho], [1, 1])\n",
    "betaV = (1 + rho*dalphardT + rho*T*d2alphardTdrho)*(rho*R)\n",
    "print(f'betaV: {betaV} Pa/K')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
